IODS course project

---

About the project

This course teaches students to use the statistical programming software R to prepare and analyze data.

Github Repository

https://github.com/blrobert/IODS-project

---

Regression and model validation

This week I have learned to setup, perform and analyze simple linear regression with multiple variables.

• Describe your work and results clearly.
• Assume the reader has an introductory course level understanding of writing and reading R code as well as statistical methods
• Assume the reader has no previous knowledge of your data or the more advanced methods you are using

lrn14<-read.csv("learning2014.csv",row.names=1)
str(lrn14)

## 'data.frame':    166 obs. of  7 variables:
##  $ gender  : Factor w/ 2 levels "F","M": 1 2 1 2 2 1 2 1 2 1 ...
##  $ Age     : int  53 55 49 53 49 38 50 37 37 42 ...
##  $ attitude: num  3.7 3.1 2.5 3.5 3.7 3.8 3.5 2.9 3.8 2.1 ...
##  $ deep    : num  3.58 2.92 3.5 3.5 3.67 ...
##  $ stra    : num  3.38 2.75 3.62 3.12 3.62 ...
##  $ surf    : num  2.58 3.17 2.25 2.25 2.83 ...
##  $ Points  : int  25 12 24 10 22 21 21 31 24 26 ...

dim(lrn14)

## [1] 166   7

This dataset includes 166 rows and 7 columns (variables). The variables are Age, Gender, (exam) Points, Attitude, Surface Learning, Strategic Learning and Deep Learning.

install.packages(“GGally”) install.packages(“ggplot2”)

library(ggplot2)

## Warning: package 'ggplot2' was built under R version 3.3.2

library(GGally)

## Warning: package 'GGally' was built under R version 3.3.2

pairs(lrn14[-1],col=lrn14$gender)

ggpairs(lrn14, mapping=aes(col=gender, alpha=0.3), lower=list(combo=wrap("facethist",bins=20)))

summary(lrn14)

##  gender       Age           attitude          deep            stra      
##  F:110   Min.   :17.00   Min.   :1.400   Min.   :1.583   Min.   :1.250  
##  M: 56   1st Qu.:21.00   1st Qu.:2.600   1st Qu.:3.333   1st Qu.:2.625  
##          Median :22.00   Median :3.200   Median :3.667   Median :3.188  
##          Mean   :25.51   Mean   :3.143   Mean   :3.680   Mean   :3.121  
##          3rd Qu.:27.00   3rd Qu.:3.700   3rd Qu.:4.083   3rd Qu.:3.625  
##          Max.   :55.00   Max.   :5.000   Max.   :4.917   Max.   :5.000  
##       surf           Points     
##  Min.   :1.583   Min.   : 7.00  
##  1st Qu.:2.417   1st Qu.:19.00  
##  Median :2.833   Median :23.00  
##  Mean   :2.787   Mean   :22.72  
##  3rd Qu.:3.167   3rd Qu.:27.75  
##  Max.   :4.333   Max.   :33.00

Most of the variables have relatively normal distributions, with the exception of age, which has a right-skewed distribution. The male attitude distribution could also be considered left-skewed. Most of the variables have weak relationships with each other, as exhibited by correlations less than 0.1. Notably, points and attitude have a strong positive relationship with a correlation value of 0.437. Deep learning and surface learning have a fairly strong negative correlation at -0.324, meaning that as deep learning increases, the tendency toward surface learning decreases. Surprisingly, median deep learning is higher than median surface learning, although there is higher variation with regard to deep learning.

my_model<-lm(Points~attitude+deep+surf, data=lrn14)
summary(my_model)

## 
## Call:
## lm(formula = Points ~ attitude + deep + surf, data = lrn14)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -16.9168  -3.1487   0.3667   3.8326  11.3519 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  18.3551     4.7124   3.895 0.000143 ***
## attitude      3.4661     0.5766   6.011 1.18e-08 ***
## deep         -0.9485     0.7903  -1.200 0.231815    
## surf         -1.0911     0.8360  -1.305 0.193669    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.313 on 162 degrees of freedom
## Multiple R-squared:  0.2024, Adjusted R-squared:  0.1876 
## F-statistic:  13.7 on 3 and 162 DF,  p-value: 5.217e-08

The only tested explanatory variable that exhibits a statistically significant relationship with the target variable points is attitude, as indicated by a p-value less than 0.05. For every 1 point increase in a student’s measured attitude, there is a 6.011 point increase in points scored on the exam. There is no significant relationship different from 0 between the explanatory variables deep learning and surface learning and the target variable points (exam score).

Ordinary least squares regression determines the Betas, or best fit lines, that minimize the differences between the fitted value and the observed value for each of the explanatory variables. The provided t-values represent the slopes of the linear best fit lines for each explanatory variable, and thus the relationship betewen the explanatory variable and the target variable. The p-value determines whether that is relationship is statistically significant, or due to chance.

my_model2<-lm(Points~attitude,data=lrn14)

summary(my_model2)

## 
## Call:
## lm(formula = Points ~ attitude, data = lrn14)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -16.9763  -3.2119   0.4339   4.1534  10.6645 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  11.6372     1.8303   6.358 1.95e-09 ***
## attitude      3.5255     0.5674   6.214 4.12e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.32 on 164 degrees of freedom
## Multiple R-squared:  0.1906, Adjusted R-squared:  0.1856 
## F-statistic: 38.61 on 1 and 164 DF,  p-value: 4.119e-09

Removing the explanatory variables surface learning and deep learning, which did not exhibit a statistically significant relationship with the target variable points, attitude remains statistically significant and exhibits a slightly larger relationship with/effect on points. The intercept also increases in value from 3.895 to 6.358. The median residual increases from the prior model, while the multiple r-squared value falls from 0.2024 to 0.1906, meaning that the second model has less explanatory power than the first.

R-squared is a measure of how close the true observations are to the fitted line, and thus a determination of how well the model explains the variation of observations around the data’s mean. It tells you how much of the variance in the target variable can be explained by the explanatory variables. R-squared values range from 0 to 1, with 0 indicating that the model explains none of the variation, and 1 indicating that the explanatory variables used explain all of the variation in the target variable. My model has an R-squared of 0.2024, meaning that it has relatively low explanatory power and there are substantial differences between the fitted lines for each explanatory variable and the actual observations.

plot(my_model2, which=c(1:2,5))

The assumptions of the linear regression model are that the errors are normally distributed, the errors are not correlated and the errors have constant varaiance (and, as such, the size of an error does not depend on the explanatory variables). For my model, the normality assumption is reasonable, because the points are for the most part tightly fitted to the line in the Q-Q Plot and there are no extreme outliers. The constant variance assumption is also reasonable, because there is a random spread of points with no clear pattern in the graph comparing residuals and model predictions. The Residuals vs. Leverage plot indicates that the data has regular leverage, because there is no significant outlier (the leverage values range from 0 to 0.04) that has a disproportionate effect on the model.

---

Logistic Regression

Introduction to logistic regression and cross validation

alc <- read.csv("alc.csv",row.names=1)

colnames(alc)

##  [1] "school"     "sex"        "age"        "address"    "famsize"   
##  [6] "Pstatus"    "Medu"       "Fedu"       "Mjob"       "Fjob"      
## [11] "reason"     "nursery"    "internet"   "guardian"   "traveltime"
## [16] "studytime"  "failures"   "schoolsup"  "famsup"     "paid"      
## [21] "activities" "higher"     "romantic"   "famrel"     "freetime"  
## [26] "goout"      "Dalc"       "Walc"       "health"     "absences"  
## [31] "G1"         "G2"         "G3"         "alc_use"    "high_use"

dim(alc)

## [1] 382  35

The data inclues 382 rows (students) and 35 variables/columns. Each variable indicates attributes and grade information about students in either math or Portuguese language courses.

I will test Pstatus, studytime, failures and higher in relation to high/low alcohol consumption. For those students with Pstatus apart (false), there will be a higher probability of high alcohol consumption. For those students with higher weekly study times, there will be a lower probability of high aclohol consumption. For those students with a higher number of failures, there will be a higher probability of high alcohol consumption. For those students who want to pursue higher education, there will be a lower probablility of high aclohol consumption.

install.packages(“GGally”) install.packages(“ggplot2”)

library(ggplot2)
library(GGally)

bar<-ggplot(data=alc, aes(x=high_use))
bar+facet_wrap("Pstatus")+geom_bar()

bar+facet_wrap("studytime")+geom_bar()

bar+facet_wrap("failures")+geom_bar()

bar+facet_wrap("higher")+geom_bar()

table(high_use=alc$high_use, Pstatus=alc$Pstatus)

##         Pstatus
## high_use   A   T
##    FALSE  26 242
##    TRUE   12 102

table(high_use=alc$high_use, studytime=alc$studytime)

##         studytime
## high_use   1   2   3   4
##    FALSE  58 135  52  23
##    TRUE   42  60   8   4

table(high_use=alc$high_use, failures=alc$failures)

##         failures
## high_use   0   1   2   3
##    FALSE 244  12  10   2
##    TRUE   90  12   9   3

table(high_use=alc$high_use, higher=alc$higher)

##         higher
## high_use  no yes
##    FALSE   9 259
##    TRUE    9 105

Based on the bar graphs and cross-tabulations, my hypotheses generally seem to be appropriate. For students whose parents live together, a higher proportion have low alcohol consumption compared with students whose parents are separated. Generally, the more students study the more likely they are to have low alcohol consumption. For students who plan to pursue higher education, a higher proportion have low alcohol consumption compared to students who don’t plan to pursue higher education. In terms of failures, however, there is a less clear relationship with alcohol consumption.

model<-glm(high_use~Pstatus+studytime+failures+higher, data=alc, family="binomial")
summary(model)

## 
## Call:
## glm(formula = high_use ~ Pstatus + studytime + failures + higher, 
##     family = "binomial", data = alc)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.5570  -0.7981  -0.7981   1.3000   2.0759  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  0.41983    0.68428   0.614 0.539520    
## PstatusT    -0.08093    0.37705  -0.215 0.830043    
## studytime   -0.52540    0.15683  -3.350 0.000808 ***
## failures     0.34848    0.19764   1.763 0.077866 .  
## higheryes   -0.26871    0.53569  -0.502 0.615936    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 465.68  on 381  degrees of freedom
## Residual deviance: 443.89  on 377  degrees of freedom
## AIC: 453.89
## 
## Number of Fisher Scoring iterations: 4

coef(model)

## (Intercept)    PstatusT   studytime    failures   higheryes 
##  0.41983005 -0.08093225 -0.52540334  0.34848048 -0.26871220

The only stastically significant predictor of high aclohol consumption at the 0 level is study time. For every one unit increase in study time, there is a 0.52540 decresase in a student’s log odds of exhibiting high alcohol consumption.

install.packages(“tidyr”)

library(tidyr)

## Warning: package 'tidyr' was built under R version 3.3.2

OR<-coef(model)%>%exp
CI<-exp(confint(model))

## Waiting for profiling to be done...

cbind(OR,CI)

##                    OR     2.5 %    97.5 %
## (Intercept) 1.5217029 0.3876426 5.7868685
## PstatusT    0.9222562 0.4488646 1.9917762
## studytime   0.5913168 0.4308370 0.7978098
## failures    1.4169129 0.9625861 2.1018616
## higheryes   0.7643632 0.2675755 2.2483346

Pstatus is associated with almost 1:1 odds of high alcohol consumption, but the odds of success are slightly lower for students whose parents are together, meaning that PstatusT is negatively associated with success(high aclohol consumption). Study time is associated with approximately 3:5 odds, indicating that those students who study more have lower odds of exhibiting high alcohol consumption. Failures is associated with approximately 3:2 odds, meaning the odds of success are higher for those students with more failures and thus failure is positively associated with high alcohol consumption. Higher education exhibits 3:4 odds, indicating that the desire to puruse higher education is negatively associated with high alchol consumption. Since none of the Confidence Intervals contain 1.0, all of the explanatory variables exhibit statistically significant associations with the target variable (high alcohol consumption) at the p<0.05 level.

Overall, the relationships posited in my hypothesis were correct, but only study time exhibited a statistically significant relationship with alcohol consumption in the regression model.

model2<-glm(high_use~studytime, data=alc, family="binomial")
summary(model2)

## 
## Call:
## glm(formula = high_use ~ studytime, family = "binomial", data = alc)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.0603  -0.8314  -0.8314   1.2993   2.1010  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   0.3209     0.3076   1.043    0.297    
## studytime    -0.6029     0.1530  -3.941  8.1e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 465.68  on 381  degrees of freedom
## Residual deviance: 448.31  on 380  degrees of freedom
## AIC: 452.31
## 
## Number of Fisher Scoring iterations: 4

library(dplyr)

## 
## Attaching package: 'dplyr'

## The following object is masked from 'package:GGally':
## 
##     nasa

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

probabilities<-predict(model2, type="response")
alc<-mutate(alc, probability=probabilities)
alc<-mutate(alc, prediction=probability>0.5)
table(high_use=alc$high_use, prediction=alc$prediction)

##         prediction
## high_use FALSE
##    FALSE   268
##    TRUE    114

select(alc, studytime, high_use, probability, prediction) %>% head(20)

##    studytime high_use probability prediction
## 1          2    FALSE   0.2921823      FALSE
## 2          2    FALSE   0.2921823      FALSE
## 3          2     TRUE   0.2921823      FALSE
## 4          3    FALSE   0.1842725      FALSE
## 5          2    FALSE   0.2921823      FALSE
## 6          2    FALSE   0.2921823      FALSE
## 7          2    FALSE   0.2921823      FALSE
## 8          2    FALSE   0.2921823      FALSE
## 9          2    FALSE   0.2921823      FALSE
## 10         2    FALSE   0.2921823      FALSE
## 11         2    FALSE   0.2921823      FALSE
## 12         3    FALSE   0.1842725      FALSE
## 13         1    FALSE   0.4299752      FALSE
## 14         2    FALSE   0.2921823      FALSE
## 15         3    FALSE   0.1842725      FALSE
## 16         1    FALSE   0.4299752      FALSE
## 17         3    FALSE   0.1842725      FALSE
## 18         2    FALSE   0.2921823      FALSE
## 19         1     TRUE   0.4299752      FALSE
## 20         1    FALSE   0.4299752      FALSE

Interestingly, using the model with only study time as an explanatory variable, there are no predicted probabilities above .5 and thus no predictions of true (high alcohol consumption).

library(ggplot2)
g<-ggplot(alc, aes(x=probability, y=high_use, col=prediction))
g+geom_point()

loss_func<-function(class,prob){
  n_wrong<-abs(class-prob)>0.5
  mean(n_wrong)
}

loss_func(class=alc$high_use, prob=alc$probability)

## [1] 0.2984293

The training error is approximately .3, meaning that 30 percent of students are inaccurately classified in terms of their alcohol consumption (high/low). Simply guessing without a strategy gives you a 50 percent chance of being correct in a binary situation such as this, and guessing with some intuition/strategy most likely gives you an even higher chance of being correct. Thus, the performance of the model is perceivably not that much different from a simple guessing strategy.

install.packages(“boot”)

library(boot)
cv<-cv.glm(data=alc, cost=loss_func, glmfit=model, K=10)
cv$delta[1]

## [1] 0.3141361

cv<-cv.glm(data=alc, cost=loss_func, glmfit=model2, K=10)
cv$delta[1]

## [1] 0.2984293

Both models have higher average prediction error than the DataCamp model. To find such a model, you could test different amounts and types of explanatory variables together.

---

Clustering and Classification

library(MASS)

## 
## Attaching package: 'MASS'

## The following object is masked from 'package:dplyr':
## 
##     select

data("Boston")
str(Boston)

## 'data.frame':    506 obs. of  14 variables:
##  $ crim   : num  0.00632 0.02731 0.02729 0.03237 0.06905 ...
##  $ zn     : num  18 0 0 0 0 0 12.5 12.5 12.5 12.5 ...
##  $ indus  : num  2.31 7.07 7.07 2.18 2.18 2.18 7.87 7.87 7.87 7.87 ...
##  $ chas   : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ nox    : num  0.538 0.469 0.469 0.458 0.458 0.458 0.524 0.524 0.524 0.524 ...
##  $ rm     : num  6.58 6.42 7.18 7 7.15 ...
##  $ age    : num  65.2 78.9 61.1 45.8 54.2 58.7 66.6 96.1 100 85.9 ...
##  $ dis    : num  4.09 4.97 4.97 6.06 6.06 ...
##  $ rad    : int  1 2 2 3 3 3 5 5 5 5 ...
##  $ tax    : num  296 242 242 222 222 222 311 311 311 311 ...
##  $ ptratio: num  15.3 17.8 17.8 18.7 18.7 18.7 15.2 15.2 15.2 15.2 ...
##  $ black  : num  397 397 393 395 397 ...
##  $ lstat  : num  4.98 9.14 4.03 2.94 5.33 ...
##  $ medv   : num  24 21.6 34.7 33.4 36.2 28.7 22.9 27.1 16.5 18.9 ...

dim(Boston)

## [1] 506  14

The dataset Boston has 506 rows and 14 columns/variables. The variables are characterised by factors that influence housing values in the suburbs of Boston.

install.packages(“corrplot”)

library(tidyr)
library(corrplot)
pairs(Boston)

cor_matrix<-cor(Boston)%>%round
cor_matrix

##         crim zn indus chas nox rm age dis rad tax ptratio black lstat medv
## crim       1  0     0    0   0  0   0   0   1   1       0     0     0    0
## zn         0  1    -1    0  -1  0  -1   1   0   0       0     0     0    0
## indus      0 -1     1    0   1  0   1  -1   1   1       0     0     1    0
## chas       0  0     0    1   0  0   0   0   0   0       0     0     0    0
## nox        0 -1     1    0   1  0   1  -1   1   1       0     0     1    0
## rm         0  0     0    0   0  1   0   0   0   0       0     0    -1    1
## age        0 -1     1    0   1  0   1  -1   0   1       0     0     1    0
## dis        0  1    -1    0  -1  0  -1   1   0  -1       0     0     0    0
## rad        1  0     1    0   1  0   0   0   1   1       0     0     0    0
## tax        1  0     1    0   1  0   1  -1   1   1       0     0     1    0
## ptratio    0  0     0    0   0  0   0   0   0   0       1     0     0   -1
## black      0  0     0    0   0  0   0   0   0   0       0     1     0    0
## lstat      0  0     1    0   1 -1   1   0   0   1       0     0     1   -1
## medv       0  0     0    0   0  1   0   0   0   0      -1     0    -1    1

corrplot(cor_matrix, method="circle",type="upper", cl.pos="b", tl.pos="d", tl.cex = 0.6 )

According to the correlation plot, there appear to be strong relationships between a number of the variables in the dataset. For example, there is a strong, positive relationship between per capita crime rates and property tax rates. There is a strong, negative relationship between property tax rates and average distances to five Boston employment centers. Both of these relationships seem counterintuitive.

boston_scaled<-scale(Boston)
summary(boston_scaled)

##       crim                 zn               indus        
##  Min.   :-0.419367   Min.   :-0.48724   Min.   :-1.5563  
##  1st Qu.:-0.410563   1st Qu.:-0.48724   1st Qu.:-0.8668  
##  Median :-0.390280   Median :-0.48724   Median :-0.2109  
##  Mean   : 0.000000   Mean   : 0.00000   Mean   : 0.0000  
##  3rd Qu.: 0.007389   3rd Qu.: 0.04872   3rd Qu.: 1.0150  
##  Max.   : 9.924110   Max.   : 3.80047   Max.   : 2.4202  
##       chas              nox                rm               age         
##  Min.   :-0.2723   Min.   :-1.4644   Min.   :-3.8764   Min.   :-2.3331  
##  1st Qu.:-0.2723   1st Qu.:-0.9121   1st Qu.:-0.5681   1st Qu.:-0.8366  
##  Median :-0.2723   Median :-0.1441   Median :-0.1084   Median : 0.3171  
##  Mean   : 0.0000   Mean   : 0.0000   Mean   : 0.0000   Mean   : 0.0000  
##  3rd Qu.:-0.2723   3rd Qu.: 0.5981   3rd Qu.: 0.4823   3rd Qu.: 0.9059  
##  Max.   : 3.6648   Max.   : 2.7296   Max.   : 3.5515   Max.   : 1.1164  
##       dis               rad               tax             ptratio       
##  Min.   :-1.2658   Min.   :-0.9819   Min.   :-1.3127   Min.   :-2.7047  
##  1st Qu.:-0.8049   1st Qu.:-0.6373   1st Qu.:-0.7668   1st Qu.:-0.4876  
##  Median :-0.2790   Median :-0.5225   Median :-0.4642   Median : 0.2746  
##  Mean   : 0.0000   Mean   : 0.0000   Mean   : 0.0000   Mean   : 0.0000  
##  3rd Qu.: 0.6617   3rd Qu.: 1.6596   3rd Qu.: 1.5294   3rd Qu.: 0.8058  
##  Max.   : 3.9566   Max.   : 1.6596   Max.   : 1.7964   Max.   : 1.6372  
##      black             lstat              medv        
##  Min.   :-3.9033   Min.   :-1.5296   Min.   :-1.9063  
##  1st Qu.: 0.2049   1st Qu.:-0.7986   1st Qu.:-0.5989  
##  Median : 0.3808   Median :-0.1811   Median :-0.1449  
##  Mean   : 0.0000   Mean   : 0.0000   Mean   : 0.0000  
##  3rd Qu.: 0.4332   3rd Qu.: 0.6024   3rd Qu.: 0.2683  
##  Max.   : 0.4406   Max.   : 3.5453   Max.   : 2.9865

summary(Boston)

##       crim                zn             indus            chas        
##  Min.   : 0.00632   Min.   :  0.00   Min.   : 0.46   Min.   :0.00000  
##  1st Qu.: 0.08204   1st Qu.:  0.00   1st Qu.: 5.19   1st Qu.:0.00000  
##  Median : 0.25651   Median :  0.00   Median : 9.69   Median :0.00000  
##  Mean   : 3.61352   Mean   : 11.36   Mean   :11.14   Mean   :0.06917  
##  3rd Qu.: 3.67708   3rd Qu.: 12.50   3rd Qu.:18.10   3rd Qu.:0.00000  
##  Max.   :88.97620   Max.   :100.00   Max.   :27.74   Max.   :1.00000  
##       nox               rm             age              dis        
##  Min.   :0.3850   Min.   :3.561   Min.   :  2.90   Min.   : 1.130  
##  1st Qu.:0.4490   1st Qu.:5.886   1st Qu.: 45.02   1st Qu.: 2.100  
##  Median :0.5380   Median :6.208   Median : 77.50   Median : 3.207  
##  Mean   :0.5547   Mean   :6.285   Mean   : 68.57   Mean   : 3.795  
##  3rd Qu.:0.6240   3rd Qu.:6.623   3rd Qu.: 94.08   3rd Qu.: 5.188  
##  Max.   :0.8710   Max.   :8.780   Max.   :100.00   Max.   :12.127  
##       rad              tax           ptratio          black       
##  Min.   : 1.000   Min.   :187.0   Min.   :12.60   Min.   :  0.32  
##  1st Qu.: 4.000   1st Qu.:279.0   1st Qu.:17.40   1st Qu.:375.38  
##  Median : 5.000   Median :330.0   Median :19.05   Median :391.44  
##  Mean   : 9.549   Mean   :408.2   Mean   :18.46   Mean   :356.67  
##  3rd Qu.:24.000   3rd Qu.:666.0   3rd Qu.:20.20   3rd Qu.:396.23  
##  Max.   :24.000   Max.   :711.0   Max.   :22.00   Max.   :396.90  
##      lstat            medv      
##  Min.   : 1.73   Min.   : 5.00  
##  1st Qu.: 6.95   1st Qu.:17.02  
##  Median :11.36   Median :21.20  
##  Mean   :12.65   Mean   :22.53  
##  3rd Qu.:16.95   3rd Qu.:25.00  
##  Max.   :37.97   Max.   :50.00

The numbers in the dataset were scaled down overall, as was the range.

boston_scaled<-as.data.frame(boston_scaled)
scaled_crim<-boston_scaled$crim
quants<-quantile(scaled_crim)
summary(quants)

##      Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
## -0.419400 -0.410600 -0.390300  1.742000  0.007389  9.924000

string<-c("low", "med_low", "med_high", "high")
crime<-cut(scaled_crim, breaks=quants, include.lowest=TRUE, label=string)
summary(crime)

##      low  med_low med_high     high 
##      127      126      126      127

boston_scaled<-dplyr::select(boston_scaled, -crim)
boston_scaled<-data.frame(boston_scaled, crime)

n<-nrow(boston_scaled)
samp<-sample(n, size=n*0.8)
train<-boston_scaled[samp,]
test<-boston_scaled[-samp,]

lda.fit<-lda(crime~.,data=train)
lda.fit

## Call:
## lda(crime ~ ., data = train)
## 
## Prior probabilities of groups:
##       low   med_low  med_high      high 
## 0.2673267 0.2400990 0.2376238 0.2549505 
## 
## Group means:
##                  zn      indus        chas        nox          rm
## low       0.8885998 -0.9011685 -0.12651054 -0.8651577  0.43361013
## med_low  -0.1566000 -0.2888065  0.01179157 -0.5406093 -0.17857441
## med_high -0.3970195  0.2184795  0.13778554  0.4106588  0.03564202
## high     -0.4872402  1.0170891 -0.08120770  1.0454118 -0.48033094
##                 age        dis        rad        tax     ptratio
## low      -0.8695139  0.8665520 -0.6894375 -0.7292938 -0.42596927
## med_low  -0.2355538  0.2879135 -0.5509001 -0.4766305 -0.05231797
## med_high  0.4638692 -0.3910414 -0.4339567 -0.3119223 -0.15315640
## high      0.8046973 -0.8630158  1.6384176  1.5142626  0.78111358
##                black       lstat        medv
## low       0.37408443 -0.76101461  0.50635596
## med_low   0.31651592 -0.09374727 -0.01701815
## med_high  0.04108174  0.07675086  0.07197759
## high     -0.71226428  0.90701275 -0.70418858
## 
## Coefficients of linear discriminants:
##                 LD1         LD2         LD3
## zn       0.08150335  0.61623334 -0.95768361
## indus    0.04537990 -0.21365597  0.20903420
## chas    -0.11761334 -0.05013481  0.11248816
## nox      0.30203491 -0.74367359 -1.56110424
## rm      -0.13821916 -0.05915772 -0.23417902
## age      0.29592046 -0.36834064 -0.03840641
## dis     -0.05669083 -0.16030236 -0.01774081
## rad      3.42129172  0.97315017 -0.08034830
## tax      0.02549263  0.01035384  0.67604958
## ptratio  0.08061452 -0.04385477 -0.41161506
## black   -0.09242340  0.06197783  0.14752305
## lstat    0.17219088 -0.15568157  0.28561737
## medv     0.14994846 -0.32084803 -0.22701803
## 
## Proportion of trace:
##    LD1    LD2    LD3 
## 0.9538 0.0342 0.0120

lda.arrows<-function(x, myscale=1, arrow_heads=0.1, color="red", tex=0.75, choices=c(1,2)){
  heads<-coef(x)
  arrows(x0=0, y0=0,
    x1=myscale*heads[,choices[1]],
    y1=myscale*heads[,choices[2]],
col=color, length=arrow_heads)
text(myscale*heads[,choices],labels=row.names(heads),
      cex=tex,col=color, pos=3)
}
classes<-as.numeric(train$crime)
plot(lda.fit, dimen=2, col=classes, pch=classes)
lda.arrows(lda.fit, myscale=1)

Low, med_low and med_high training data are all intermixed and spacially tight, while the high training data are mostly clustered together, separate from other categories and spacially distant (with the exception of a few med_high data intermixed and one high entry mixed into the other cluster). It appears as though the high predictions will be the most accurate (most likely because the training data are the most spacially distant from the others), followed by the low predictions. The med_low and med_high predictions seem as though they will be relatively inaccurate, and the categories could most likely be slimmed to two or three classifications/clusters.

correct_classes<-test$crime
test<-dplyr::select(test,-crime)
lda.pred<-predict(lda.fit, newdata=test)
table(correct=correct_classes, predicted=lda.pred$class)

##           predicted
## correct    low med_low med_high high
##   low       14       4        1    0
##   med_low   10      16        3    0
##   med_high   1       8       19    2
##   high       0       0        0   24

The LDA model was much more accurate at predicting the med_high and high results-it correctly predicted every high result. Low results were predicted correctly just over 50 percent of the time, although the incorrect predictions were all med_low. Med_low results were predicted under 50 percent of the time, and predictions ranged form low to med_high, though med_low was predicted the most.

library(MASS)
data("Boston")
boston_scaled<-scale(Boston)
dist_eu<-dist(boston_scaled)
summary(dist_eu)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.1343  3.4620  4.8240  4.9110  6.1860 14.4000

km<-kmeans(dist_eu, centers=15)
pairs(Boston, col=km$cluster)

k_max<-15
twcss<-sapply(1:k_max, function(k){kmeans(dist_eu,k)$tot.withinss})
plot(1:k_max, twcss, type='b')

km2<-kmeans(dist_eu, centers=2)
pairs(Boston, col=km2$cluster)

Average number of rooms per dwelling (rm), nitrogen oxide concentration (nox) and proportion of owner-occupied units built prior to 1940 (age) seem to affect the clustering results.

library(MASS)
data("Boston")
boston_scaled<-scale(Boston)
dist_eu<-dist(boston_scaled)
km3<-kmeans(dist_eu, centers=4)

boston_scaled<-as.data.frame(boston_scaled)
lda.fit2<-lda(km3$cluster~.,data=boston_scaled)
lda.fit2

## Call:
## lda(km3$cluster ~ ., data = boston_scaled)
## 
## Prior probabilities of groups:
##         1         2         3         4 
## 0.2272727 0.4308300 0.2114625 0.1304348 
## 
## Group means:
##         crim         zn      indus       chas        nox         rm
## 1  0.2797949 -0.4872402  1.1892663 -0.2723291  0.8998296 -0.2770011
## 2 -0.3894453 -0.2173896 -0.5212959 -0.2723291 -0.5203495 -0.1157814
## 3 -0.3912182  1.2671159 -0.8754697  0.5739635 -0.7359091  0.9938426
## 4  1.4330759 -0.4872402  1.0689719  0.4435073  1.3439101 -0.7461469
##          age        dis        rad        tax     ptratio       black
## 1  0.7716696 -0.7723199  0.9006160  1.0311612  0.60093343 -0.01717546
## 2 -0.3256000  0.3182404 -0.5741127 -0.6240070  0.02986213  0.34248644
## 3 -0.6949417  0.7751031 -0.5965444 -0.6369476 -0.96586616  0.34190729
## 4  0.8575386 -0.9620552  1.2941816  1.2970210  0.42015742 -1.65562038
##        lstat        medv
## 1  0.6116223 -0.54636549
## 2 -0.2813666 -0.01314324
## 3 -0.8200275  1.11919598
## 4  1.1930953 -0.81904111
## 
## Coefficients of linear discriminants:
##                 LD1        LD2         LD3
## crim    -0.18113078 -0.5012256 -0.60535205
## zn      -0.43297497 -1.0486194  0.67406151
## indus   -1.37753200  0.3016928  1.07034034
## chas     0.04307937 -0.7598229 -0.22448239
## nox     -1.04674638 -0.3861005 -0.33268952
## rm       0.14912869 -0.1510367  0.67942589
## age      0.09897424  0.0523110  0.26285587
## dis     -0.13139210 -0.1593367 -0.03487882
## rad     -0.65824136  0.5189795  0.48145070
## tax     -0.28903561 -0.5773959  0.10350513
## ptratio -0.22236843  0.1668597 -0.09181715
## black    0.42730704  0.5843973  0.89869354
## lstat   -0.24320629 -0.6197780 -0.01119242
## medv    -0.21961575 -0.9485829 -0.17065360
## 
## Proportion of trace:
##    LD1    LD2    LD3 
## 0.7596 0.1768 0.0636

lda.arrows<-function(x, myscale=1, arrow_heads=0.1, color="red", tex=0.75, choices=c(1,2)){
  heads<-coef(x)
  arrows(x0=0, y0=0,
    x1=myscale*heads[,choices[1]],
    y1=myscale*heads[,choices[2]],
col=color, length=arrow_heads)
text(myscale*heads[,choices],labels=row.names(heads),
      cex=tex,col=color, pos=3)
}
classes=as.numeric(km3$cluster)
plot(lda.fit2, dimen=2, col=classes, pch=classes)
lda.arrows(lda.fit2, myscale=3)

Property tax rates (tax), proportion of residential land zoned (zn) and nitrogen oxides concentration (nox) are the most influential lienar separators for the clusters. These variables mostly influence the fourth cluster.

---

Dimensionality Reduction Techniques

human<-human<-read.csv("create_human.csv",row.names=1)
str(human)

## 'data.frame':    155 obs. of  8 variables:
##  $ edu_ratio  : num  97 94.5 95.8 96 89.1 ...
##  $ labor_ratio: num  65 65.3 68.3 62.5 64.5 ...
##  $ edu_exp    : num  17.5 20.2 15.8 18.7 17.9 16.5 18.6 16.5 15.9 19.2 ...
##  $ life_exp   : num  81.6 82.4 83 80.2 81.6 80.9 80.9 79.1 82 81.8 ...
##  $ GNI        : int  64992 42261 56431 44025 45435 43919 39568 52947 42155 32689 ...
##  $ MMR        : int  4 6 6 5 6 7 9 28 11 8 ...
##  $ ABR        : num  7.8 12.1 1.9 5.1 6.2 3.8 8.2 31 14.5 25.3 ...
##  $ fem_rep    : num  39.6 30.5 28.5 38 36.9 36.9 19.9 19.4 28.2 31.4 ...

dim(human)

## [1] 155   8

Human is a dataset with 155 rows/countries and 8 columns/variables. The variables are key measures in the Human Development and Gender Inequality Indices. edu_ratio (mutated from data in the Gender Inequality Index) describes the ratio of female and male populations with secondary education in each country. labor_ratio (mutated from the GII data) describes ratio of labor force participation of females and males in each country. edu_exp (from the Human Development Index) describes the expected years of schooling of the general population. life_exp (from the HDI) describes the life expectancy of the general population. GNI (HDI) describes the Gross National Income per capita of a country in 2014. MMR (GII) describes the maternal mortality ratio (deaths per 100,000 live births), and ABR (GII) describes the adolescent birth rate (births per 1,000 women ages 15-19). fem_rep (GII) describes the percentage of seats in parliament held by women.

install.packages(“tidyr”) install.packages(“ggplot2”) install.packages(“GGally”) install.packages(“corrplot”)

library(tidyr)
library(ggplot2)
library(GGally)
library(corrplot)
ggpairs(human)

cor(human)%>%corrplot

Life expectancy and adolescent birth rate have a relatively strong negative correlation, meaning that as the adolescent birth rate increases, life expectancy decreases, and vice versa. Expected education and education ratio have a relatively strong positive correlation, indicating that the ratio of female to male populations with secondary education increases and thus female educational inequality decreases as expected education increases, and vice versa. GNI, MMR and ABR all have right-skewed distriubtions, while labor_ratio, edu_exp and fem_rep has relatively normal distributions. Edu_ratio and life_exp have left-skewed distributions.

summary(human)

##    edu_ratio       labor_ratio       edu_exp         life_exp    
##  Min.   :  2.05   Min.   :40.90   Min.   : 5.40   Min.   :49.00  
##  1st Qu.: 33.27   1st Qu.:56.65   1st Qu.:11.25   1st Qu.:66.30  
##  Median : 56.95   Median :63.40   Median :13.50   Median :74.20  
##  Mean   : 58.03   Mean   :63.45   Mean   :13.18   Mean   :71.65  
##  3rd Qu.: 86.65   3rd Qu.:69.30   3rd Qu.:15.20   3rd Qu.:77.25  
##  Max.   :100.00   Max.   :89.15   Max.   :20.20   Max.   :83.50  
##       GNI              MMR              ABR            fem_rep     
##  Min.   :   581   Min.   :   1.0   Min.   :  0.60   Min.   : 0.00  
##  1st Qu.:  4198   1st Qu.:  11.5   1st Qu.: 12.65   1st Qu.:12.40  
##  Median : 12040   Median :  49.0   Median : 33.60   Median :19.30  
##  Mean   : 17628   Mean   : 149.1   Mean   : 47.16   Mean   :20.91  
##  3rd Qu.: 24512   3rd Qu.: 190.0   3rd Qu.: 71.95   3rd Qu.:27.95  
##  Max.   :123124   Max.   :1100.0   Max.   :204.80   Max.   :57.50

pca_human<-prcomp(human)
s<-summary(pca_human)
pca_pr<-round(100*s$importance[2,],digits=1)
pc_lab<-paste0(names(pca_pr),"(",pca_pr,"%)")
biplot(pca_human, cex=c(0.8,1), col=c("grey40","deeppink2"), xlab=pc_lab[1],ylab=pc_lab[2])

## Warning in arrows(0, 0, y[, 1L] * 0.8, y[, 2L] * 0.8, col = col[2L], length
## = arrow.len): zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, y[, 1L] * 0.8, y[, 2L] * 0.8, col = col[2L], length
## = arrow.len): zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, y[, 1L] * 0.8, y[, 2L] * 0.8, col = col[2L], length
## = arrow.len): zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, y[, 1L] * 0.8, y[, 2L] * 0.8, col = col[2L], length
## = arrow.len): zero-length arrow is of indeterminate angle and so skipped

In the nonstandardized PCA, there appears to be no clear correlation between the original features themselves and the features and the principal components. The only exception is GNI, which exhibits a strong correlation with the first principal component. Interestingly, the first principal component explains 100 percent of the variation in the data, possibly because GNI explains all the variation. Most of the countries appear to be clustered in the middle, with little variation, explaining why the original features seem to have little relation to the principal component and each other and why PC1 explains all of the variability in the data.

human_std<-scale(human)
summary(human_std)

##    edu_ratio         labor_ratio           edu_exp       
##  Min.   :-1.91047   Min.   :-2.200017   Min.   :-2.7378  
##  1st Qu.:-0.84475   1st Qu.:-0.663486   1st Qu.:-0.6782  
##  Median :-0.03672   Median :-0.004972   Median : 0.1140  
##  Mean   : 0.00000   Mean   : 0.000000   Mean   : 0.0000  
##  3rd Qu.: 0.97695   3rd Qu.: 0.570617   3rd Qu.: 0.7126  
##  Max.   : 1.43259   Max.   : 2.507135   Max.   : 2.4730  
##     life_exp            GNI               MMR               ABR         
##  Min.   :-2.7188   Min.   :-0.9193   Min.   :-0.6992   Min.   :-1.1325  
##  1st Qu.:-0.6425   1st Qu.:-0.7243   1st Qu.:-0.6496   1st Qu.:-0.8394  
##  Median : 0.3056   Median :-0.3013   Median :-0.4726   Median :-0.3298  
##  Mean   : 0.0000   Mean   : 0.0000   Mean   : 0.0000   Mean   : 0.0000  
##  3rd Qu.: 0.6717   3rd Qu.: 0.3712   3rd Qu.: 0.1932   3rd Qu.: 0.6030  
##  Max.   : 1.4218   Max.   : 5.6890   Max.   : 4.4899   Max.   : 3.8344  
##     fem_rep       
##  Min.   :-1.8203  
##  1st Qu.:-0.7409  
##  Median :-0.1403  
##  Mean   : 0.0000  
##  3rd Qu.: 0.6127  
##  Max.   : 3.1850

pca_human_std<-prcomp(human_std)

s2<-summary(pca_human_std)
pca_pr2<-round(100*s2$importance[2,],digits=1)
pc_lab2<-paste0(names(pca_pr2),"(",pca_pr2,"%)")
biplot(pca_human_std, cex=c(0.8,1), col=c("grey40","deeppink2"), xlab=pc_lab2[1],ylab=pc_lab2[2])

summary(pca_human_std)

## Importance of components:
##                           PC1    PC2     PC3     PC4     PC5     PC6
## Standard deviation     2.1585 1.0679 0.88106 0.68579 0.59186 0.53987
## Proportion of Variance 0.5824 0.1426 0.09703 0.05879 0.04379 0.03643
## Cumulative Proportion  0.5824 0.7249 0.82195 0.88074 0.92453 0.96096
##                            PC7     PC8
## Standard deviation     0.44297 0.34071
## Proportion of Variance 0.02453 0.01451
## Cumulative Proportion  0.98549 1.00000

In the standardized PCA, the correlations are more evident and the countries are more spread out with greater variation. Scaling the data down allows for greater visualization of variability in the data.

Female representation in parliament and the ratio of labor force participation betewen females and males are relatively correlated with each other and with principal component 2 (fem_rep has the stronger correlation with PC2), indicating that PC2 most likely explains the variation in the country data relating to gender inequality. However, PC2 only captures 14 percent of the variability in the data.The rest of the variables are strongly correlated with principal component one and with each other, although in different directions. Maternal Mortality Ratio and Adolescent Birth Rate are strongly correlated with each other and PC1. Gross National Income per capita, the ratio of female and male populations with secondary education in each country, expected years of education, and life expectancy are strongly correlated with each other and PC1, but in the opposite direction of MMR and ABR. It is likely that these variables explain the human (or economic) development variability between the countries expressed in the data. High ABR and MMR contribute negatively to human development while higher GNI, education expectations, life expectancies and levels of female secondary educations point to more economically developed countries. Labor_ratio seems to be equally correlated with PC1 and PC2, which makes sense because higher female participation in the labor force indicates both lower gender inequality and higher levels of economic development. PC1 explains 58.2 percent of the variability in the data.

Overall, PC1 and PC2 indicate the influence of different types of original features on the overall data. However, PC1 has significantly more explanatory power and is highly correlated with more original features. The features in PC1 are more related to human/economic development, and thus Principal Component One describes the variation in the data based on human/economic development. PC2 is more related to variation between countries in terms of gender inequality, and has less ability to capture the variation in the data based on the reality that more of the original features are related to the HDI. Even variables from the GII have a relationship with PC1, because development is related to gender equality.

install.packages(“FactoMineR”) install.packages(“dplyr”)

library(FactoMineR)

## Warning: package 'FactoMineR' was built under R version 3.3.2

data("tea")
str(tea)

## 'data.frame':    300 obs. of  36 variables:
##  $ breakfast       : Factor w/ 2 levels "breakfast","Not.breakfast": 1 1 2 2 1 2 1 2 1 1 ...
##  $ tea.time        : Factor w/ 2 levels "Not.tea time",..: 1 1 2 1 1 1 2 2 2 1 ...
##  $ evening         : Factor w/ 2 levels "evening","Not.evening": 2 2 1 2 1 2 2 1 2 1 ...
##  $ lunch           : Factor w/ 2 levels "lunch","Not.lunch": 2 2 2 2 2 2 2 2 2 2 ...
##  $ dinner          : Factor w/ 2 levels "dinner","Not.dinner": 2 2 1 1 2 1 2 2 2 2 ...
##  $ always          : Factor w/ 2 levels "always","Not.always": 2 2 2 2 1 2 2 2 2 2 ...
##  $ home            : Factor w/ 2 levels "home","Not.home": 1 1 1 1 1 1 1 1 1 1 ...
##  $ work            : Factor w/ 2 levels "Not.work","work": 1 1 2 1 1 1 1 1 1 1 ...
##  $ tearoom         : Factor w/ 2 levels "Not.tearoom",..: 1 1 1 1 1 1 1 1 1 2 ...
##  $ friends         : Factor w/ 2 levels "friends","Not.friends": 2 2 1 2 2 2 1 2 2 2 ...
##  $ resto           : Factor w/ 2 levels "Not.resto","resto": 1 1 2 1 1 1 1 1 1 1 ...
##  $ pub             : Factor w/ 2 levels "Not.pub","pub": 1 1 1 1 1 1 1 1 1 1 ...
##  $ Tea             : Factor w/ 3 levels "black","Earl Grey",..: 1 1 2 2 2 2 2 1 2 1 ...
##  $ How             : Factor w/ 4 levels "alone","lemon",..: 1 3 1 1 1 1 1 3 3 1 ...
##  $ sugar           : Factor w/ 2 levels "No.sugar","sugar": 2 1 1 2 1 1 1 1 1 1 ...
##  $ how             : Factor w/ 3 levels "tea bag","tea bag+unpackaged",..: 1 1 1 1 1 1 1 1 2 2 ...
##  $ where           : Factor w/ 3 levels "chain store",..: 1 1 1 1 1 1 1 1 2 2 ...
##  $ price           : Factor w/ 6 levels "p_branded","p_cheap",..: 4 6 6 6 6 3 6 6 5 5 ...
##  $ age             : int  39 45 47 23 48 21 37 36 40 37 ...
##  $ sex             : Factor w/ 2 levels "F","M": 2 1 1 2 2 2 2 1 2 2 ...
##  $ SPC             : Factor w/ 7 levels "employee","middle",..: 2 2 4 6 1 6 5 2 5 5 ...
##  $ Sport           : Factor w/ 2 levels "Not.sportsman",..: 2 2 2 1 2 2 2 2 2 1 ...
##  $ age_Q           : Factor w/ 5 levels "15-24","25-34",..: 3 4 4 1 4 1 3 3 3 3 ...
##  $ frequency       : Factor w/ 4 levels "1/day","1 to 2/week",..: 1 1 3 1 3 1 4 2 3 3 ...
##  $ escape.exoticism: Factor w/ 2 levels "escape-exoticism",..: 2 1 2 1 1 2 2 2 2 2 ...
##  $ spirituality    : Factor w/ 2 levels "Not.spirituality",..: 1 1 1 2 2 1 1 1 1 1 ...
##  $ healthy         : Factor w/ 2 levels "healthy","Not.healthy": 1 1 1 1 2 1 1 1 2 1 ...
##  $ diuretic        : Factor w/ 2 levels "diuretic","Not.diuretic": 2 1 1 2 1 2 2 2 2 1 ...
##  $ friendliness    : Factor w/ 2 levels "friendliness",..: 2 2 1 2 1 2 2 1 2 1 ...
##  $ iron.absorption : Factor w/ 2 levels "iron absorption",..: 2 2 2 2 2 2 2 2 2 2 ...
##  $ feminine        : Factor w/ 2 levels "feminine","Not.feminine": 2 2 2 2 2 2 2 1 2 2 ...
##  $ sophisticated   : Factor w/ 2 levels "Not.sophisticated",..: 1 1 1 2 1 1 1 2 2 1 ...
##  $ slimming        : Factor w/ 2 levels "No.slimming",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ exciting        : Factor w/ 2 levels "exciting","No.exciting": 2 1 2 2 2 2 2 2 2 2 ...
##  $ relaxing        : Factor w/ 2 levels "No.relaxing",..: 1 1 2 2 2 2 2 2 2 2 ...
##  $ effect.on.health: Factor w/ 2 levels "effect on health",..: 2 2 2 2 2 2 2 2 2 2 ...

dim(tea)

## [1] 300  36

library(dplyr)
keep_columns<-c("Tea", "How", "how", "sugar", "where", "lunch")
tea_time<-dplyr::select(tea, one_of(keep_columns))
str(tea_time)

## 'data.frame':    300 obs. of  6 variables:
##  $ Tea  : Factor w/ 3 levels "black","Earl Grey",..: 1 1 2 2 2 2 2 1 2 1 ...
##  $ How  : Factor w/ 4 levels "alone","lemon",..: 1 3 1 1 1 1 1 3 3 1 ...
##  $ how  : Factor w/ 3 levels "tea bag","tea bag+unpackaged",..: 1 1 1 1 1 1 1 1 2 2 ...
##  $ sugar: Factor w/ 2 levels "No.sugar","sugar": 2 1 1 2 1 1 1 1 1 1 ...
##  $ where: Factor w/ 3 levels "chain store",..: 1 1 1 1 1 1 1 1 2 2 ...
##  $ lunch: Factor w/ 2 levels "lunch","Not.lunch": 2 2 2 2 2 2 2 2 2 2 ...

gather(tea_time)%>%ggplot(aes(value))+geom_bar()+facet_wrap("key", scales="free")+theme(axis.text.x=element_text(angle=45,hjust=1,size=8))

## Warning: attributes are not identical across measure variables; they will
## be dropped

mca<-MCA(tea_time, graph=FALSE)
summary(mca)

## 
## Call:
## MCA(X = tea_time, graph = FALSE) 
## 
## 
## Eigenvalues
##                        Dim.1   Dim.2   Dim.3   Dim.4   Dim.5   Dim.6
## Variance               0.279   0.261   0.219   0.189   0.177   0.156
## % of var.             15.238  14.232  11.964  10.333   9.667   8.519
## Cumulative % of var.  15.238  29.471  41.435  51.768  61.434  69.953
##                        Dim.7   Dim.8   Dim.9  Dim.10  Dim.11
## Variance               0.144   0.141   0.117   0.087   0.062
## % of var.              7.841   7.705   6.392   4.724   3.385
## Cumulative % of var.  77.794  85.500  91.891  96.615 100.000
## 
## Individuals (the 10 first)
##                       Dim.1    ctr   cos2    Dim.2    ctr   cos2    Dim.3
## 1                  | -0.298  0.106  0.086 | -0.328  0.137  0.105 | -0.327
## 2                  | -0.237  0.067  0.036 | -0.136  0.024  0.012 | -0.695
## 3                  | -0.369  0.162  0.231 | -0.300  0.115  0.153 | -0.202
## 4                  | -0.530  0.335  0.460 | -0.318  0.129  0.166 |  0.211
## 5                  | -0.369  0.162  0.231 | -0.300  0.115  0.153 | -0.202
## 6                  | -0.369  0.162  0.231 | -0.300  0.115  0.153 | -0.202
## 7                  | -0.369  0.162  0.231 | -0.300  0.115  0.153 | -0.202
## 8                  | -0.237  0.067  0.036 | -0.136  0.024  0.012 | -0.695
## 9                  |  0.143  0.024  0.012 |  0.871  0.969  0.435 | -0.067
## 10                 |  0.476  0.271  0.140 |  0.687  0.604  0.291 | -0.650
##                       ctr   cos2  
## 1                   0.163  0.104 |
## 2                   0.735  0.314 |
## 3                   0.062  0.069 |
## 4                   0.068  0.073 |
## 5                   0.062  0.069 |
## 6                   0.062  0.069 |
## 7                   0.062  0.069 |
## 8                   0.735  0.314 |
## 9                   0.007  0.003 |
## 10                  0.643  0.261 |
## 
## Categories (the 10 first)
##                        Dim.1     ctr    cos2  v.test     Dim.2     ctr
## black              |   0.473   3.288   0.073   4.677 |   0.094   0.139
## Earl Grey          |  -0.264   2.680   0.126  -6.137 |   0.123   0.626
## green              |   0.486   1.547   0.029   2.952 |  -0.933   6.111
## alone              |  -0.018   0.012   0.001  -0.418 |  -0.262   2.841
## lemon              |   0.669   2.938   0.055   4.068 |   0.531   1.979
## milk               |  -0.337   1.420   0.030  -3.002 |   0.272   0.990
## other              |   0.288   0.148   0.003   0.876 |   1.820   6.347
## tea bag            |  -0.608  12.499   0.483 -12.023 |  -0.351   4.459
## tea bag+unpackaged |   0.350   2.289   0.056   4.088 |   1.024  20.968
## unpackaged         |   1.958  27.432   0.523  12.499 |  -1.015   7.898
##                       cos2  v.test     Dim.3     ctr    cos2  v.test  
## black                0.003   0.929 |  -1.081  21.888   0.382 -10.692 |
## Earl Grey            0.027   2.867 |   0.433   9.160   0.338  10.053 |
## green                0.107  -5.669 |  -0.108   0.098   0.001  -0.659 |
## alone                0.127  -6.164 |  -0.113   0.627   0.024  -2.655 |
## lemon                0.035   3.226 |   1.329  14.771   0.218   8.081 |
## milk                 0.020   2.422 |   0.013   0.003   0.000   0.116 |
## other                0.102   5.534 |  -2.524  14.526   0.197  -7.676 |
## tea bag              0.161  -6.941 |  -0.065   0.183   0.006  -1.287 |
## tea bag+unpackaged   0.478  11.956 |   0.019   0.009   0.000   0.226 |
## unpackaged           0.141  -6.482 |   0.257   0.602   0.009   1.640 |
## 
## Categorical variables (eta2)
##                      Dim.1 Dim.2 Dim.3  
## Tea                | 0.126 0.108 0.410 |
## How                | 0.076 0.190 0.394 |
## how                | 0.708 0.522 0.010 |
## sugar              | 0.065 0.001 0.336 |
## where              | 0.702 0.681 0.055 |
## lunch              | 0.000 0.064 0.111 |

plot(mca, invisible=c("ind"), habillage="quali")

summary(tea$where)

##          chain store chain store+tea shop             tea shop 
##                  192                   78                   30

summary(tea$how)

##            tea bag tea bag+unpackaged         unpackaged 
##                170                 94                 36

According to the MCA summary, how and where contribute substantially to dim 1 and moreso than the other variables to dimension 2.

Dim 1 explains 15 percent of the variance in the data, while Dim 2 explains 14 percent. Many of the categories are similar, with the exception of other, chain store+tea shop, tea bag+unpackaged, green, unpackaged and tea shop. Most of these categories are from the variables how and where, indicating their influence on the first two dimensions and on the variation in the data as a whole.

---